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Abstract. The nonaxisymmctric Ekman flow excited inside a neutron star following 
a rotational glitch is calculated analytically including stratification and compressibility. 
For the largest glitches, the gravitational wave strain produced by the hydrodynamic 
£/3 . mass quadrupole moment approaches the sensitivity range of advanced long-baseline 

| interferometers. It is shown that the viscosity, compressibility, and orientation of the 

■ star can be inferred in principle from the width and amplitude ratios of the Fourier 

peaks (at the spin frequency and its first harmonic) observed in the gravitational 
' wave spectrum in the + and x polarizations. These transport coefficients constrain 

H ' the equation of state of bulk nuclear matter, because they depend sensitively on the 

Ch ' degree of superfluidity 

PACS numbers: 04.30.Db, 95.30.Lz, 95.30.Sf, 97.60. Gb, 97.60.Jd 

<n : 

m . , 

fT) . 1. Introduction 

^ : 

Glitches are tiny, impulsive angular accelerations observed in rotation-powered pulsars. 
In total, 287 glitches have been discovered in 101 radio pulsars over four decades of radio 
pulse timing experiments, the majority since the advent of the Parkes Multibeam Survey, 
[1] improved multifrequency ephemerides, and better interference rejection algorithms. 
[2, 3, 4] Glitches also occur occasionally in anomalous X-ray pulsars, which are believed 
to be ultramagnetized neutron stars (magnetars). [5, 6] The glitch population has been 
analyzed statistically from several standpoints. [7, 8, 9, 10, 2, 3, 4, 11, 12] The fractional 
increase in angular velocity spans a wide range, 10" 11 < 5tt/tt < 10" 4 , suggesting that 
the angular momentum stored in the superfiuid core of the neutron star is transferred 
erratically to the solid crust (~ 1% of the total moment of inertia) via a series of 
discrete coupling events of varying strength. [13, 12] While the trigger is unknown, it is 
often ascribed to collective unpinning of quantized superfiuid vortices, e.g. when crust- 
superfluid differential rotation (and hence the Magnus force) exceeds a threshold, [14] 
or when a starquake occurs, heating the crust and boosting the rate of vortex creep, [15] 
although there are dissenting views. [16] Vortex unpinning is followed by an exponential 
relaxation phase, during which crust-superfluid corotation is restored by viscous forces. 
[17, 18] 

It is natural to speculate that glitches excite nonaxisymmetric motions, generating 
gravitational radiation. There are several plausible scenarios in which this can happen. 
(1) A bulk two-stream instability occurs in the core, mediated by entrainment between 
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the neutron superfluid and charged species. [19] (2) A surface two-stream instability 
occurs at the interface between the 1 So- and 3 P 2 -paired neutron superfluids. [20] (3) 
The crust cracks, tilts, and precesses. [21, 11] (4) Crust-core differential rotation 
drives meridional circulation and, at high Reynolds number, superfluid turbulence. 
[18, 22, 23, 11, 24] (5) Vortex unpinning excites pulsation modes in the mult i- component 
core superfluid, e.g. acoustic, gravity, and Rossby waves, [25, 26] and quasiradial 
oscillations. [27, 28] On the basis of general arguments, invoking the conservation 
of energy and angular momentum, the gravitational wave strain is predicted to reach 
h ~ 1CT 26 from pulsars that glitch strongly and frequently, like Vela. [26, 27] 

In this paper, we calculate analytically the gravitational wave signal generated 
during the relaxation phase of a glitch, when the fluid interior of a neutron star responds 
to an impulsive, nonaxisymmetric angular acceleration of the crust. We idealize the 
spin-up event as an enduring step increase in the angular velocity of a semi-rigid, 
fluid-filled container, deliberately simplifying the initial and boundary conditions to 
keep the focus on how effects like stratification and compressibility (see below) show 
up in gravitational-wave observables. (Besides, the detailed microphysics of the spin- 
up event is unknown in most of the scenarios listed in the previous paragraph.) For 
example, we treat the interior as a single Newtonian fluid with an average viscosity 
and density, whereas in reality it is a multi-component superfluid with viscous and 
inviscid components, whose properties vary with depth. We therefore neglect the two- 
stream dynamics in scenarios (1) and (2), nor do we allow for an interface between the 
inner crust and outer core [cf. scenario (2)]. The circulation of a rotating superfluid 
is carried by quantized vortices, which can pin metastably in between nuclear lattice 
sites, or at defects (e.g., dislocations, grain boundaries), in the crust. [29, 30] Yet we 
neglect the effects of vortex pinning in our model [cf. scenario (5)], except to note that 
pinning introduces an ambiguity in the boundary conditions, which is not yet resolved 
in experiments with terrestrial superfluids, let alone in a neutron star. (We expand on 
this issue in §2.4 below.) Finally, we do not allow the crust to precess [cf. scenario 
(3)], partly to keep the analysis tractable, and partly because precession of the coupled 
crust-core system is tied up intimately with the strength of vortex pinning [31], which 
we neglect. (In general, the crust possesses a misaligned mass quadrupole moment, 
whose gravitational wave signal must be added independently to the glitch signal we 
compute.) Regrettably, there is no published literature on the nonlinear evolution of 
the instabilities which underlie the above glitch scenarios, whether multi-fluid modes in 
(l)-(4) or vortex Kelvin waves in (5); the relevant numerical simulations have not yet 
been performed. Yet the nonlinear saturation state of the instability determines the 
degree of nonaxisymmetry (and hence the gravitational wave spectrum) of the spun up 
container in reality. [32, 33] This shortcoming is felt most keenly in scenario (4), where 
numerical simulations are known to reveal nonaxisymmetric meridional circulation even 
in the small-shear limit [18, 22], an effect excluded artificially by the linear analysis in 
this paper. 

With these approximations in mind, we construct a toy spin-up model in §2 by 
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solving the nonaxisymmetric Ekman problem using the traditional method of multiple 
scales, generalizing previous axisymmetric analyses. [34, 35] We next investigate the 
effects of stratification (§3) and compressibility (§4). In §5, we compute the gravitational 
wave strain generated by the mass quadrupole. The current quadrupole, which is 
typically smaller, is evaluated in a forthcoming paper. In §6, §7, we show that the 
amplitude, frequency, phase, quality factor, and polarization can be inverted to measure 
the viscosity and compressibility of the stellar interior and hence infer its state of 
superfluidity. Gravitational radiation from glitches therefore probes directly the bulk 
properties of quantum condensates at nuclear densities ~ 10 15 gcur 3 . The latter topic is 
attracting growing attention, both theoretically [36, 37, 38, 39, 40, 41] and in heavy-ion 
collision experiments [42, 43, 44, 45]. 

2. Nonaxisymmetric Ekman flow after a glitch 

During a glitch, the angular velocity of the stellar crust increases by a small amount 
80, < fi, on a time-scale that is unresolved by observations [46] (cf. some Crab glitches 
[47]). In reality, the angular acceleration of the crust is likely to be nonaxisymmetric 
to some degree, as argued in Section 1. Consequently, the resulting spin-up (Ekman) 
flow is likely to be nonaxisymmetric as well. In this section, we calculate the Ekman 
flow analytically, generalizing previous axisymmetric analyses. [48, 35] The exact flow 
pattern is affected strongly by stratification and compressibility. These properties, which 
directly probe the physics of bulk nuclear matter, therefore leave their imprint on the 
resulting gravitational wave signal. 

The initial and boundary conditions adopted to describe the spin-up event 
hydrodynamically are discussed critically in §2.4 and §2.6. We adopt take the standard 
form favored in the neutron star literature [34, 35, 18, 49]: a permanent, discontinuous 
spin-up of a fluid-filled container, which is semi-rigid to accommodate nonaxisymmetry. 
These conditions have a long history of being used successfully as the simplest way 
to capture the spirit of a glitch in a toy model. Nevertheless, they are plainly not 
self-consistent in three important respects, discussed furthur in §2.4 and §2.6. (1) The 
physics of the glitch trigger is unknown, but it is likely to involve vortex pinning, which 
these hydrodynamic conditions neglect. (2) The superfluid interior couples frictionally 
to the solid crust in a way that allows partial slippage at the interface [50, 51, 52]. The 
degree of slippage is still not characterised fully in laboratory experiments with liquid 
helium, let alone in neutron stars. (3) The back reaction from vortex turbulence and 
magnetic fluxoids on the Ekman boundary layer, which contributes to the boundary 
conditions on the interior flow, remains an unsolved problem [18, 23, 24, 33, 53, 54, 32]. 
In view of these unavoidable uncertainties, the idealizations in the model (e.g. semi- 
rigidity) are tolerable. 

For simplicity, we treat the spin up of a cylinder, following Ref. [35]. The 
nonaxisymmetric Ekman problem in a sphere has not yet been solved analytically in 
general, although some progress has been made in special limits. [55] 
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2.1. Model equations 

Consider a cylinder of height 2L and radius L, containing a compressible Newtonian 
fluid with uniform kinematic viscosity v. (Note that the effective (eddy averaged) v 
may be boosted substantially by turbulent Reynolds stresses, as in terrestrial oceanic 
or atmospheric flows [56, 57]). Initially, the cylinder rotates with angular velocity 
Q = Qe z about the z axis. In the rotating frame, the compressible Navier-Stokes 
equation, including Coriolis and centrifugal forces, takes the form 

^+v-Vv + 2tixv = --Vp + g + vV 2 v+^V(V -v)+V (ln 2 r 2 ) .(1) 
at p 3 V2 / 

The symbols v, p, p, and g denote the fluid velocity, fluid density, pressure, and 

gravitational acceleration respectively. Magnetic forces are neglected for the sake of 

simplicity, even though they are known to modify the spin up of the core (which 

is coupled to the crust by cyclotron- vortex waves). [34, 58] We adopt cylindrical 

coordinates (r, 0, z) in the rotating frame. The gravitational acceleration is taken to 

be uniform and perpendicular to the midplane of the cylinder, viz. 

f -ge z if^>0 

y \ +ge z if*<0, 1 1 

where g is constant. This form of the gravitational field is standard in both the neutron 
star literature [34, 35] and the classic papers of geomechanics [59, 48, 55]. It is sourced 
by a singular surface mass distribution at z = and is therefore nominally unphysical. 
However, it leads to a self-consistent spin-up model, renders the problem analytically 
tractable, and has been validated by comparisons with numerical simulations involving 
smooth mass distributions [60]. Henceforth, we restrict attention to the region z > 0, 
as the flow is symmetric about the midplane of the cylinder. 
The fluid satisfies the continuity equation 

^ + V-(pv) = 0. (3) 

The energy equation is written in a form that relates the convective derivatives of the 
pressure and density, viz. 

4g+-v)p, (4) 

where c 2 is the propagation speed of an acoustic disturbance in the Lagrangian frame, 
determined by the equation of state (e.g., adiabatic or isothermal). 

It is useful to recast (2)-(4) in dimensionless form, when preparing to separate the 
slow and fast dynamics of the flow. The characteristic time-scale for Ekman pumping 
is t E = E' 1 / 2 ^}- 1 , where we define the dimensionless Ekman number 

The characteristic length-scale is L. The velocity scale for the spin-up flow is L5Q. The 
remaining scales are chosen so that the hydrostatic forces dominate inertial forces in the 
momentum equation. The pressure scale is taken to be p c = pogL, the density scale p 
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is taken to be the equilibrium density in the midplane z — 0. These choices lead to the 
following scaled equations: 



eF [E 1/2 ^ + ev-Vv + 2e z xv 



= --Vp-e z + eFE 
P 



V 2 v + \v(V-v) 



fFV(V) , (G) 



El/2 ft + eV "(H = o, (7) 

(E l ' 2 j t + ev ■ vjp = K (V/ 2 | + ev-Vjp. (8) 

Equations (6)-(8) introduce three dimensionless quantities: the Rossby number e = 
5Q/Q, the Froude Number F = LQ 2 /g, and the scaled compressibility K = gL/c 2 . 

At time t — 0, the cylinder is accelerated impulsively. If the Rossby number is small, 
we can solve for the equilibrium and spin-up (Ekman) flows separately by making the 
perturbation expansion p i— > p + e5p, p i— > p + e5p. We also replace v with 5v since it is 
of order e. In the rotating frame, the equilibrium flow is steady and axisymmetric, with 
p = p(r, z), p = p(r, z) and 

= --\/p-e z + FV (-r 2 ) . (9) 



P V2 

On the other hand, the spin-up flow is unsteady and nonaxisymmetric, with Sp 

5p(r, 0, z, t),6p = 8p(r, 0, z, t),5v = 8v(r, 0, z, t), 



= —-V5p — —e z + FE 
P P 



V 2 5v + -V (V • Sv) 



and 



+ ^V(^), (10) 

^ 1/2 ^ + V-( P ^) = 0, (11) 

E 1 ^ + 6v-V P = k(e 1 ' 2 ?!£ + W-Vp\ . (12) 
at \ at J 

The effects of stratification and compressibility enter the spin-up flow through (12), 

which is sometimes written alternatively as [61] 

E^ = KE^- P 5v.A. (13) 

dt dt ' v 1 

The Schwarzchild discriminant, 

7 l Vp + Ke z -KFv( l -r 2 ] , (14) 



p r V2 
characterises the buoyancy. The buoyant restoring force experienced by a fluid element 
when displaced from equilibrium in the z direction is proportional to —gd(ln p)/dz — 
g 2 /c 2 . It increases as the density gradient steepens and decreases as the compressibility 
increases. The frequency of oscillation of the fluid element is called the Brunt- Vaisala 
frequency. When scaled to Q, the Brunt- Vaisala frequency is related to the Schwarzchild 

1 /2 

discriminant by N = {—A z /F) 
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2. 2. Equilibrium flow 

To solve (9) and (10)-(12) analytically, we follow Refs [35] and [48] and neglect 
centrifugal terms, as the Froude number is typically small in a neutron star. However, 
we retain all other terms containing F until the end of the calculation, as FN 2 and K 
can be of similar magnitude. We also assume that p~ x dpjdz is uniform throughout the 
star to simplify the problem analytically. The stratification of the equilibrium flow can 
arise from either compositional or entropic gradients. 
Solving (9) we obtain 



p(z) 



-K 3 z 



K7 1 e~ KsZ 



(15) 
(16) 



where the dimensionless quantity, K s = L/z s = —Lp~ 1 dp/dz is proportional to the 
reciprocal of the stratification length scale, z s . The density at the midplane (z = 0) is 
chosen to be the fudicial density scale p , as discussed in §2.1. The pressure at the surface 
[z = L) can be made to vanish by adding a constant to (16), but we leave it nonzero 
here (without loss of generality) in order to maintain consistency with the assumption of 
uniform sound speed, e.g. if the flow is adiabatic, we have c 2 = dp/ dp = •yp/p = ^gz s , 
where 7 is the ratio of specific heats. 



2.3. Spin-up flow 

Neglecting the centrifugal force terms, the vector components of (10)-(12) are 
F ( E i/2dSv r 



dt 



25v, 



X 



r 2 J r 2 d(p 3 dr 



(V-^17,) 



F E 1 / 



^ dSv *+25v r )=- 1 -^( S A+FE 



dt 



x 



o<p \ p J 

' 2 1 \ . 2 d5v r 1 d _ . J n , 



and 



FE 1 ' 2 



dSv z 
~dt 



x 



P P 



1 d 



V^ + --(V.^) 



E 1 ' 2 



(19) 
(20) 

(21) 
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where N 2 = (K s — K) / F. We now proceed according to the method of multiple 
scales by expanding Sv, 5p and Sp in perturbation series in the small parameter E x l 2 . 
[62, 48, 59, 35] To order O(E ), we find 

(22) 

(23) 

(24) 
(25) 

8z ' (26) 
where the superscript '0' denotes a zeroth-order quantity. In view of (22)-(24), the 
0(E°) continuity equation (25) is identically satisfied. To order 0(E 1 ^ 2 ), we find 

i J_d® J_d_ (5_p^\ 

Vr AF 8r 2Frd<P\ p J ' 1 ] 

5v * = AFr^ + 2Fd-r{—) ' (28) 
1 <9$ 

| (^) + V ■ Sv' = KM , (30) 

where we define $ = —d(5p°/p)/dt, and the superscript '1' denotes a first-order 
quantity. Upon substituting (22)-(29) into the 0(E 1 ^ 2 ) continuity equation (30) we 
arrive at 

1 d ( <9$\ 1 <9 2 $ 4K s d<S> 4 <9 2 $ 

+ 1^1^ = 0. (31) 



r dr \ dr ) r 2 d(p 2 N 2 dz N 2 dz 2 

Equation (31) is solvable by separation of variables. It is slightly more general than the 
equation for $ in Ref. [35] as it allows for the possibility that FN 2 and K are of similar 
magnitude, which is likely in a neutron star. 



2.4. Boundary conditions 

The boundary conditions on $ are determined by the boundary conditions on the 
velocity fields, given by equations (22) and (23). Note that the boundary conditions are 
imposed on the 0(E°) flow. To impose boundary conditions on the 0(E 1 ^ 2 ) flow, one 

must specify Sp 1 and utilise (27) and (28). The general solution to (31) has the form 

00 00 

$(r, (j), z, t) = F J m {\ mn r) [A mn (t) cos(m0) + B mn (t) sin(m^)] Z mn (z) ,(32) 

m=0 n=l 

where m > is an integer, X mn is the n-th root of J m (X) = 0, and A mn and B mn 
are coefficients which are functions of time. The prefactor F is extracted explicitly in 
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anticipation that $ is a quantity of order F. The solution is periodic in 0, regular in 
the limit r — > 0, and vanishes at the cylinder wall r = 1. 

To specify the boundary conditions on the axial component of the velocity, we 
require that the flow be symmetric about the midplane of the cylinder, viz. 5v\ = at 
z — 0. The final boundary condition on Z mn is arbitrary, and we choose Z mn (l) = 1. 
Hence we find 



The boundary conditions for the axisymmetric mode (m = 0) at r = 1 differ from 
those of previous authors [48, 35]. For the axisymmetric mode, the zeroth-order radial 
flow is zero everywhere by virtue of (22), and Refs. [48, 35] impose no penetration on 
the first-order radial flow (27), giving d&/dr — at r = 1. However, for the higher 
order nonaxisymmetric modes (m > 0), we are forced, by (22), to demand <9<£>/<90 = 
at r = 1, so that the no-penetration boundary condition on the zeroth order radial flow 
is satisfied. For generality, we impose the same boundary condition on the axisymmetric 
mode as on the nonaxisymmetric modes, namely <9<£>/<90 = at r = 1. The resultant 
solution qualitatively resembles that found in Refs. [48, 35]. 

The above boundary conditions, based on a single Newtonian fluid with an average 
viscosity, are grossly idealized. In fact, the liquid interior of the star is a multi-component 
mixture of charged, viscous fluid and uncharged, inviscid superfluid. The boundary 
conditions satisfied by the superfluid mixture are uncertain, not just in neutron stars, 
but also in terrestrial experiments with liquid helium. A large body of numerical and 
experimental research [50, 63, 52] suggests that the behaviour of the superfluid lies 
somewhere between no slip and perfect slip where it touches the container, with the 
exact amount of slip determined by the density of the (turbulent) quantized vortices 
in the boundary layer, which is poorly known. In a neutron star, this uncertainty is 
magnified by the presence of superconducting magnetized fluxoids, which modify the 
structure of the boundary layer in a nonaxisymmetric way. [33, 53, 32] Moreover, the 
velocity components at the container are partly determined by pinning, which this 
paper neglects, but which is predicted to vary in strength across macroscopic capacitive 
domains [32, 64]. 

2.5. Temporal evolution 

In Ekman pumping, the temporal evolution of the interior flow is controlled by the 
boundary conditions at the top and bottom faces of the cylinder. The mass flux into 
(out of) the thin viscous boundary layer is related to the circulation just outside the 
layer by [48, 59, 35] 




{FN 2 - p_) e^ z - (FN 2 - (3+) e?- z 
(FN 2 - /?_) e?+ - {FN 2 - (3 + ) 1F- 



(33) 



with 




(34) 



8v 



z\z=±l ~ 



T l -E l / 2 {Vx5v) z 



z=±l 



(35) 
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This relation follows from the structure of the Ekman layer in an incompressible flow, 
but it applies equally when the flow is compressible; the fourth term on the right hand 
side of (1) introduces corrections of 0(E) to (35). Upon differentiating (35) with respect 
to time and substituting (22), (23), and (32), we obtain 



dSvl 



dt 



z=±l 



= T i AL<f(r,0,±l,t) • 



(36) 



Combining (29) and (36), we find that the (m, n)-th Ekman mode relaxes exponentially 
on the Ekman time-scale as $ oc exp(— u mn t), with [35] 

(EN 2 - (3_) e?+ - (FN 2 - (3+) e?- 



\2 



^FK + \ 2 mn )(e^-eP-) 
The general solution for the pressure perturbation can now be written down: 

5p°(r,(j),z,t) Sjf^fafaz) 



(37) 



P(z) 

oo oo 

- F J m(^ mn r) (A mn cos mcft + B mn sin mcft) 



m=0 n=l 

x Z mn (z)e~ L 



lt . (38) 

In (38), the function Sp^r, (ft, z) corresponds to the steady-state pressure profile of the 
spun up cylinder. All that remains is to determine the coefficients A mn and B mn , which 
are now constants, by specifying the initial and final conditions of the spin up. 



2.6. Initial conditions 

If the pressure perturbation is zero initially, then A mn and B mn are determined by 
evaluating (38) at z — 1. Specifying the steady-state solution Sp^/p is equivalent to 
specifying the steady-state velocity fields obeying (22) and (23). We therefore require a 
steady-state, nonaxisymmetric solution to the Navier-Stokes equations, namely 

M^lil = F £ r m {r 2 _ 1} cog(m0) (39) 
Pi 1 ) ^0 

which ensures that there is no penetration at r = 1, and that the flow vanishes smoothly 
on the rotation axis. For m — 0, the familiar steady-state, axisymmetric flow is recovered 
[35, 48]. The corresponding velocity boundary conditions on the top and bottom faces 
of the cylinder are (in dimensional form) 



Sv r 



\L,t) = 



m5QL 



m—l 



5v<p(r,<ft,L,t) = 



5QL 



m—l 



(m + 2) 



smmc 



m 



cosmcft 



(40) 



(41) 



The boundary conditions on the side walls are given by (40) and (41) evaluated at r = L. 
The above choice is somewhat artificial, because it implies that the top and bottom faces 
of the cylinder do not rotate rigidly, as one might expect for the crust of a neutron star. 
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However, this shortcoming is not serious when compared with the other uncertainties 
identified at the start of §2 and in §2.4 (e.g. the degree of boundary layer slippage, and 
inhomogeneous pinning [32]). Indeed, the semi-rigid flow boundary conditions described 
by (40) and (41) is a reasonable model for the outer edge of the interior flow which is 
isolated from the rigid container by the turbulent Ekman boundary layer. The model 
captures qualitatively the main elements of the astrophysical problem: a step increase 
in angular velocity, which in turn excites a nonaxisymmetric internal flow. 
The Fourier coefficients A mn in (43)-(47) can be computed from 



the first few are A u = -0.706, A 12 = 0.154, A 21 = -0.521, and A 22 = 0.148. Given 
(39), we have B mn = 0. 

The steady-state solution (39)-(41) implies that the cylinder spins up to a 
nonaxisymmetric steady state (which continually emits gravitational radiation; see §5). 
In reality, it is indeed likely that some nonaxisymmetry survives between glitches, 
as a result of inhomogeneous vortex pinning in macroscopic domains. [54, 47, 32] 
As an alternative, however, one can consider a scenario in which nonaxisymmetric 
modes are excited by the glitch trigger at t — (as discussed in §1) and subsequently 
dissipate as t — > oo, leaving only the axisymmetric mode. This scenario can be treated 
mathematically by specifying the conditions corresponding to (38) at t — and t — > oo. 
It is deferred to future work. 

2. 7. Density and velocity fields 

With the constants of the Fourier expansion determined, the density, pressure and 
velocity can be written out explicitly using (22)-(30). Transforming out of the rotating 
frame into the observer's frame, and restoring dimensional variables, the complete 
expressions read 




(42) 



v r (r,(f),z,t) 



L 2 5Vt 



mA mn J m (X mn r/L) sinm(0 - Qt) 



oo oo 



2r 



m=0 n=l 



X 



(FN 2 - e^ z l L - (FN 2 - (3+) e?- z / L 
(FN 2 - JT) eP+ - (FN 2 - (3+) e?- 




(43) 



v<i>(r,<f>,z,t) 



fir + 



L 2 5VL ~ ~ . d 



Y 12 A ranTT [J m (A mn r/L)] COSm(0 



nt) 



m=0 n=l 



X 



(FN 2 - /?_) e p + z l L - (FN 2 - ftjQ eP- z ' L 
(FN 2 - JZ) e?+ - (FN 2 - 13+) ef- 




(44) 
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. .2.^ A rnn A mn J m {A mn r / L) 



m=0 n=l 



x cosm(0 — Ot) 



— E^/^UJmrit 



(45) 



p(r, 0, z, i) = p e + po 



XI A mn J m (\ mn r / L) cos[m(0 - fit)] 



x 



9 m =0n=l 

(FN 2 - p_e~P- z / L - (FN 2 - p + ) fi + e-P+ z ' L 



(FN 2 - P-) e?+ - (FN 2 - (3+) e?- 
x (l - e^ 1 ' 2 ^ 1 ) , 



(46) 



oo oo 



p(r, 0, z, t) = p gz s e z/Zs + p Q L 2 Vt5Vt ^ ^ A mn J m (\ rnn r/ L) cos m(<f> - Qt) 

m=0 n=l 

'(FN 2 - (3_) e -P- z ' L - (FN 2 - /3+) e~?+ z / L ~ 



x 



X 



(1 



(FN 2 - e?+ - (FN 2 - p + ) 



')■ 



(47) 



3. Stratification 



Stratification acts to restrict the penetration depth of the secondary 0(E 1 ^ 2 ) flow, 
preventing the interior of the container from spinning up completely. Because the spun 
up volume decreases, the Ekman time decreases. When a fluid element is displaced 
vertically from its equilibrium position in the hydrostatic density gradient dp/dz, it 
experiences a buoyancy force, proportional to N 2 , as discussed in §2.1. The strength 
of the buoyancy force is determined in part by the fluid element's ability to adjust its 
density to match the surroundings, which is related to its compressibility (see §4 below). 
However, for a fixed compressibility, the z-dependent factors in (43) and (44) indicate 
that the rotating fluid on the midplane lags behind the final velocity at z — ±1 as 
t — > oo. The lag increases as dp/dz increases. The interior of the cylinder eventually 
spins up when vorticity diffuses out of the boundary layer on the much longer time-scale 

En- 1 . 

To illustrate the physics of stratification, consider the special case of an 
incompressible Ekman flow. To investigate this limit, we take K — > and hence 
N 2 -> K s /F. As the Froude number is small we have K s < N 2 , f3± « ±X mn N/2, 
and hence 



\ N A N 
coth ■ 



n 2 — 2 ■ < 48) 

Equation (48) shows that, in the incompressible limit, an increase in the density gradient 
reduces the spin up time. This effect, familiar from axisymmetric flows [48], is enhanced 
for higher mode numbers m. 

Figure 1 summarizes the above physics pictorially. The four panels display contour 
plots of — Qr at t = E^^cu^l in the meridional plane — Vlt = 0, for modes m — 1,5 
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Figure 1. Steady-state azimuthal velocity contours in the meridional plane <p — VLt = 
in the rotating frame, illustrating the effects of stratification and mode number. 
Contour values, top panels: v^—Q.r = 0.8, 0.6, 0.4, 0.2, (solid lines),— 0.2, —0.4 (dashed 
lines). Bottom panels: - fir = 0.8,0.6,0.4,0.2,0.1,0.05,0 (solid lines) -0.05, -0.1 
(dashed lines). Stratification and mode number: (K s ,m) = (1, 1) (top left), (5, 1) (top 
right), (1,5) (bottom left), (5,5) (bottom right). We set F = 1 and K = 0.5 in all 
panels. 



(top, bottom), stratification K s — 1,5 (left, right), and fixed compressibility K = 0.5. 
In the right-hand panels, where dp/dz is larger, the flow on the midplane lags the flow 
at z = ±1 by more than in the left-hand panels, where dp/dz is smaller. The effect is 
more pronounced in the bottom panels, where m is larger than in the top panels. 

Figure 2 illustrates how Ekman pumping proceeds with time when five modes 
(1 < m < 5) are included with equal weighting in (39). The left-hand panels are 
for K s — 1 and correspond to snapshots at tE 1 / 2 uju = 0.1,0.3,1.0 running down the 
page. The right-hand panels show the same thing for K s = 5. Again, as time passes, the 
Ekman layer spreads from the top and sides into the interior of the cylinder, but it does 
so less rapidly when the hydrostatic density gradient is greater. The time for v</, — Qr to 
reach e -1 of its maximum value at z — ±1 and r = in the left and right-hand panels 
is t — 0.53 and 0.50 respectively (in dimensionless units). 
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Figure 2. Temporal evolution of azimuthal velocity contours in the meridional plane 
<f> — Ot = in the rotating frame, showing the effect of stratification. Five modes 
(1 < m < 5), equally weighted, are included in the solution. Contour values, left panels: 
v<j> - Or = 0.8,0.6,0.4,0.2,0 (solid lines), -0.2, -0.4(dashed lines). Right panels: 
v^-Qr = 1.0,0.8,0.6,0.4,0.2,0.1,0.05,0 (solid lines), -0.05, -0.1, -0.2, -0.4(dashed 
lines). Each column contains three snapshots, at tE^-^cjn = 0.1 (top), 0.3 (middle), 
and 1.0 (bottom). Stratification parameter: K s = 1 (left), 5 (right). We set F = 1 and 
K = 0.5 in all panels. 



4. Compressibility 

Compressibility acts in two ways. (1) It reduces the buoyancy force experienced by 
an element of fluid displaced from its equilibrium position in the hydrostatic density 
gradient. It therefore competes against the density gradient to increase both the spin- 
up time and the volume of fluid spun up. (2) It also affects the spin-up time through 
mass conservation. As the axial secondary flow moves from the top (or bottom) face 
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to the midplane through the density gradient, it is compressed, decelerating the radial 
secondary flow and increasing the spin-up time. Therefore compressibility increases the 
spin-up time even when there is no buoyancy force. This is demonstrated in the zero 
buoyancy limit N — > 0, where we have K s = K, [3 + = K s , j3_ = 0, and hence 



v r (r,(j),z,t) 



/ j. \ m—l 



X 



St 



m=0 



msinm(0 — fit) , 



(49) 



v,(r,<P,z,t)=rn + ^(l-e- El/ ^) 



x 



m=0 ^ 

; (r,0,z,t) = E^LSnl 6 



m—l 



(m + 2) ^— ^ — m 



cosm(0 — fit) , (50) 



gz/c 2 



\ e gL/c 2 _ i 



m=0 



x ^2 [t ) ( m + ■*■) cosm(0 - nt) , 



(51) 



p(r,(f),z,t) = p Q e 



gz/c 2 p L 2 Q5Q c _ gz/c 2 _ e _ E i/2 wt 



oo 
m=0 



+ 

r 
I 



cosm(0 — fit) 



p(r, 0, ^, t) = p c 2 e"^ /c2 + p L 2 fi5fie^ 2 / c2 (1 - e' sl/2 ^) 



X 



E 

m=0 



cosm 



fit) . 



(52) 



(53) 



Note that v r and are independent of z in this limit. The reciprocal of the spin-up 
time is [35] 

uj K 

a = < 54 » 

When there is no buoyancy, the spin-up time increases exponentially with the 
compressibility. All mode numbers m are equally affected. In the limit c — > oo, the 
nonaxisymmetric generalisation of Ref. [62] is obtained. 

Figure 3 summarizes the above physics pictorially. The four panels display contour 
plots of ity — fir at t = E^^lu^i in the meridional plane — fit = 0, for modes m — 1, 5 
(top, bottom), compressibility K — 2,4 (left, right), and fixed stratification K s = 5. 
In the left-hand panels, where c 2 is larger, the flow on the midplane lags the flow at 
z = ±1 by more than in the right-hand panels, where c 2 is smaller. The effect is more 
pronounced on the bottom panels, where m is larger than in the top panels. 

Figure 4 illustrates how Ekman pumping proceeds with time when five modes 
(1 < m < 5) are included with equal weighting in (39). The left-hand panels are 
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Figure 3. Steady-state azimuthal velocity contours in the meridional plane 
— ilt = in the rotating frame, illustrating the effects of compressibility and 
mode number. Contour values, top panels: — fir = 0.5,0.2,0.1,0.05,0.01,0 
(solid lines), —0.01, —0.05, — O.f, —0.2 (dashed lines). Bottom panels: — fir = 
0.5,0.2,0.1,0.05,0.01,0.001,0 (solid lines), -0.001, -0.01, -0.05, -O.f , -0.2 (dashed 
lines). Compressibility and mode number: (K,m) = (2, f) (top left), (4, f) (top right), 
(2, 5) (bottom left), (4, 5) (bottom right). We set F = 1 and K = 0.5 in all panels. 



for K = 2 and correspond to snapshots at tE 1 ^ 2 Un = 0.1,0.3,1.0 running down the 
page. The right-hand panels show the same thing for K = 4. Again, as time passes, the 
Ekman layer spreads from the top and sides into the interior of the cylinder, but it does 
so more slowly when the compressibility is greater. The time for — Qr to reach e _1 
of its maximum value at z — ±1 and r = in the left and right-hand panels is t — 0.55 
and 0.63 respectively (in dimensionless units). 

5. Gravitational wave signal 

The nonaxisymmetric secondary (spin-up) flow creates a stress energy distribution with 
a time- varying mass quadrupole moment, which emits gravitational radiation. In this 
section, we calculate the gravitational wave signal from the flow solution derived in 
§2, for the special cases of polar and equatorial observers. These special cases suffice 
to demonstrate, in principle, how the signal can be inverted to infer the internal 
properties of the star, specifically its viscosity and compressibility, as described in §6. A 
more general treatment, valid for arbitrary inclination angles and including the current 
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Figure 4. Temporal evolution of azimuthal velocity contours in the meridional plane 
<fi — fit = in the rotating frame, showing the effect of compressibility. Five modes 
(1 < to < 5), equally weighted, arc included in the solution. Contour values, left 
panels: - Vlr = 0.8,0.5,0.2,0.1,0.05,0.01,0 (solid lines), -0.01,-0.05,-0.1,-0.2 
(dashed lines). Right panels: ty, - Vlr = 0.8,0.5,0.2,0.1,0.05,0.01,0.001,0 (solid 
lines), —0.001, —0.01, —0.05, —0.1, —0.2 (dashed lines). Each column contains three 
snapshots, at tE 1 l 2 LO\\ = 0.1 (top), 0.3 (middle), and 1.0 (bottom). Compressibility 
parameter: K = 2 (left), 4 (right) We set F = 1 and K s = 5 in all panels. 



quadrupole moment, is postponed to a future paper. 



5.1. Stress- energy tensor 

The relativistic stress-energy tensor of a perfect fluid is given by 



P 



(55) 
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where w M is the fluid 4-velocity and g^ v is the metric tensor. We neglect viscous terms in 
(55) as they are of order 0(E). In the quadrupole moment formalism, the gravitational 
wave signal is proportional to T 00 . Expanding T 00 in powers of v/c, expanding v, p and 
p in powers of e as in §2, and noting that we are now working in the inertial frame of 
the observer, we obtain 

T 00 = pc 2 + t[c 2 5p + 2pv-5v + 2^v-5v^ . (56) 

upto and including terms of order 0(e) and 0(v/c). Terms in (56) containing only 
equilibrium quantities are static and axisymmetric and do not contribute to the 
gravitational wave signal. The remaining terms are in the ratios 

^:2,|9.W|:2£l»-«l = l:(^):(f)' (57) 

Therefore in the limit where the centrifugal force is much smaller than the Coriolis force, 
the leading order terms are 

T 00 = 5pc 2 + 2pv-5v + 2^v- Sv (58) 

and T 00 is of characteristic magnitude QSQLpoc 2 /g. If gL/c 2 is also small, as in a 
neutron star, then the dominant stress-energy contribution is Spc 2 , with Sp given by 
(46). Note that, from equation (55) onwards, c now denotes the speed of light, not the 
sound speed (cf. §2.1). 



5.2. Wave strain for a polar observer 

In Einstein's quadrupole moment formalism, the components of the gravitational wave 
strain in the transverse traceless gauge for an observer positioned along the z' axis are 

hZl (t) = - ti™, (t) = ^ [4 V (t) - i yW (t)} (59) 

h T / y ,(t) = ^/xvW (6°) 

where an overdot denotes a time derivative, d is the distance to the source, and is 
the reduced quadrupole moment-of-inertia tensor defined by 



(61) 



Consider a polar observer located along the rotation (z) axis. Equations (59)-(61) 
can be combined with (46) to give the following wave strains in the plus and cross 
polarisations respectively: 



W*) = ?|E *2nZ 2n { [(Ail 2 - Eu 2 2n ) e~ El/2 ^ 



n 2 



4fi 2 



cos 2Vtt 



n=l 



-AE l/2 w 2n ne- El/2w ^ t sin2fit} , 



(62) 
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h(p) x (t) 



ho_ 

n 2 



E «2n^2n {4£ 1/2 U> 2n fie" £l/ W COS 2fit 
ra=l 

+ [(Ail 2 - Ecu 2 2n )e- El/2w2nt - Ail 2 ] sin2ftt} . 



(63) 



In (62) and (63), we define the dimensionless characteristic amplitude 

and the coefficients a mn = A mn J m+ i(\ mn ) / \ mn , the first few of which are an = —0.074, 
ai2 = —0.007, a 21 = —0.035, a 22 = 0.005. The factors Z mn are also dimensionless, with 



Z mn — 2 



(FN 2 - (3_)e p + - (FN 2 - (3+)e 



(FN 2 - (3J)(3-e~P- z - (FN 2 - (3 + )(3 + e- p+z 



v 2-m 



• (65) 



x / dz 

In §6, we present a simple recipe for relating features of the gravitational wave 
spectrum to the physical properties of the flow. We do this for the Fourier spectrum 
(instead of the raw signal in the time domain) to facilitate presentation of the physics; 
in practice, a time-domain template would be used to analyse real gravitational wave 
data. Noting that a 2n diminish quickly with increasing n, we can approximate (62) and 
(63) by the first terms in their series. The moduli of the Fourier transformed cross and 
plus wave strains are then 



\h 



(P)H 



MP 



+ (16uj 4 + E 2 ul l )(Eu 2 1 + u 2 



h 2 a 2 21 \z 21 \yw 



\h 



(P)x 



MP 



(AQ 2 + Eu 2 21 - uj 2 ) 2 + (2uE l l 2 u 2l ) 2 
2VL [16w 4 + 4(2fi 2 + u 2 )Eu 2 2l + E 2 ^} 



h\(x 2 21 \Z 21 \ 



| 2 /fi 4 " (AVt 2 + Euj 21 — uj 2 ) 2 + (2ujE 1 / 2 uj 2 i) 2 ' (67) 
Equations (66) and (67) describe Lorentzian spectra, as for a damped harmonic oscillator 
with resonances at oo 2 = AVt 2 + Eoo 21 and damping coefficient E 1 I 2 lo 2 i. In neutron stars, 
which have E x I 2 uj 2 \ <C 2Vl, the Fourier peaks are much narrower than their separation, 
and hence we can safely ignore the interference (product) terms between them: 

2h a 2 i\Z 21 \ 



\h(p) + (uj ~ 2Q) | = \h(p) x (uj ~ 2Q)\ 



(68) 



^(2^ - uj) 2 + Euj\ 

The spectra of the two polarizations are therefore indistinguishable for a polar observer. 

For small values of E, the approximate result (68) agrees with the exact Fourier 
transform of (62) and (63) (with terms 1 < n < 50 included) to an accuracy of A%. 
A comparison is displayed in Figure 5, where (68) is plotted together with the Fourier 
transform of (62) and (63) for E = 10~ 4 and K, N — 1. 



5. 3. Wave strain for an equatorial observer 

For an equatorial observer located along the x-axis, the plus and cross wave strains as 
functions of time are given by 



V)M) = ^ E ^Z 2n { [(An 2 - Eul) e~ El/2 ^ - Art 



2tt 2 



n=l 



cos 2Qt 



-AE 1/2 w 2n VLe 



-E 1 / 2 W2nt 



sin2fit} , 



(69) 
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Figure 5. Full plus and cross polarization spectra for a polar observer, compared with 
the approximate result for E — 10~ 4 and K, N = 1. Dashed curve: Fourier transform 
of (62); dashed-dotted curve: Fourier transform of (62); solid curve: approximate result 
(68). The sum is taken to 50 terms in (62) and (63). 



V)x(0 = 7^ E amZm {{2E 1 / 2 w ln {le~ El,2w ^ t cosfit 



n 2 



n=l 



+ [(Q 2 - Eu{ n ) e - El/2w ^ - Q 2 ] sin fit} 



(70) 



From (62) and (69), it can be seen that the plus polarisation for an equatorial observer 
is exactly half as strong as the plus polarisation for a polar observer (and n radians out 
of phase). In addition, the cross polarisation for an equatorial observer oscillates at Q 
rather than 2Q. Proceeding as in §5.2, we obtain the following (n = 1) approximations 
for the plus and cross spectra: 

h a 21 \Z 21 \ 



|^(E)xM| 



yftjii - lj) 2 + Eu 2 u ' 



(71) 



(72) 



6. Measuring the viscosity and compressibility of bulk nuclear matter 

The detailed wave form computed in §5 can be inverted to measure N, E, and K s (and 
hence K, given F. 

To illustrate how, let us consider first a polar observer. For such an observer, the 
gravitational wave spectrum peaks at u = 2Q. The amplitude of the peak is 



M 

and its width is 



h 



(p)+ 



h 



(p)+ 



M 



h 



h 



(P)x 



(P)x 



2h a 2 i\Z 21 \E- 1/2 LO 21 1 



(73) 



(74) 



A more useful expression than (73) is the product of the peak amplitude and width 



M 



h 



(P)+ 



h 



(P)+ 



= M 



h 



(P)x 



h 



(P)x 



= 2\ / 3h a 2 i\Z 21 \, 



(75) 
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which is independent of E. If the distance to the source is somehow known, along with 
Sfl (e.g. from radio timing data), then ho is known. If, in addition we can develop 
a reliable estimate of v and hence E, (74) and (75) can be solved for 1V21 and |^2i|- 
As these quantities are functions of K s and N, (74) and (75) constrain the internal 
properties of the neutron star. If ho and E are unknown, then (74) and (75) place 
constraints on the latter two quantites, as well as K s and N. 




0.1 1.0 10.0 0.1 1.0 10.0 



Figure 6. Width and amplitudes for the gravitational wave spectrum seen by 
a polar observer, showing the effect of buoyancy and stratification. Two sets of 
contours are plotted in each panel. Left: width (74), (solid curve), width x 
amplitude (75), (dashed curve). Contour values: TE' 1 / 2 = 1, 2, 5, 10, 20, 40 and 
TM/ho = 0.025,0.05,0.075,0.125,0.15,0.175,0.2. Right: width (74), (solid curve), 
amplitude (73), (dashed). Contour values: YE- 1 ' 2 = 1,2,5,10,20,40 and T/h = 
0.005, 0.025, 0.05, 0.1, 0.25, 0.5, 1. 



Figure 6 presents two panels. On the left, contours of constant V 



h 



(P)+ 



E- 1 ' 2 and 



constant Y 



h 



(P)- 



M 



(P)H 



/ho are plotted as solid and dashed curves as a function 



of buoyancy and stratification. On the right, contours of constant Y 



h 



(P)- 



E' 1 ' 2 and 



constant M /i(p)+ /h Q are plotted. The figure illustrates how measurement of the width 
and amplitude of the gravitational wave spectrum for a polar observer can be used to 
constrain the quantities YE' 1 / 2 and YM/h . In addition, if E and h are known, the 
intersection of the measured solid and dashed contours represents a direct measurement 
of the properties of bulk nuclear matter. The width of the spectral peak depends most 
on the buoyancy (nearly horizontal contours). In contrast, the amplitude of the spectral 
peak exhibits a saddle at (K s , N) ~ (1, 0.4), corresponding to an inversion in sign of the 



strain. For high buoyancies, Y 



(P)+ 



M 



h 



(P)+ 



depends only on stratification, 
e to infer the internal properties of the 



Although Figure 6 can be used in principl 
neutron star from the width and amplitude of the gravitational wave spectrum, it is 
rare in practice that h and E are known. For the special case of a polar observer, the 
plus and cross polarisations are identical, and all the gravitational radiation is emitted 
at iv — 2Q, so we only have two measurable parameters, the width and amplitude of the 
spectral peak at iv = 2Q, to constrain the unknowns K, N, ho and E. As we show below, 
an equatorial observer receives radiation at iv = Q and iv = 2Q and thereby has four 
measurable parameters at his disposal. However, an observer at an inclination angle 
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< i < 7r/2 observes signals at Q and 2fl with h + ^ h x at both frequencies. Therefore, 
such an observer has eight independent peak amplitudes and widths to play with - 
enough, indeed, to solve uniquely for K, N, ho and E, as well as the inclination angle, 
i. This general case, which is more complicated, will be investigated in a forthcoming 
paper [65]. 

For an equatorial observer, the h + spectrum peaks at 2f2, while the h x spectrum 
peaks at Q. The respective widths and amplitudes are 



^(E) + 
h(E)x 



V3E 1/2 oo 21 



V3E 1 / 2 



M 


h (E)+ 


r 




= V3h a 2 i\Z 2 i\ 


M 


h(E)x 


r 


h(E)x 


= V3h an\Z u \ 



(76) 
(77) 
(78) 
(79) 

Equations (76)-(79) contain four measureable quantities, which can be used to constrain 
the unknown quantities K, N, h and E. 

Figure 7 plots the relationships in (76)-(79) graphically. The top-left, top-right 
and bottom-left panels show contours of MT/h and T for both the cross and plus 
polarizations for different slices of parameter space. Slices are taken at E — 1CT 7 , 
N = 1, and K s — 1 respectively. These panels can be used to solve for the four unknown 
quantities K s , N, E and ho by locating the intersection of the contours corresponding 
to the measurements of T and M. The bottom-right panel has contours of the ratios 
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Because the width and 



and T 

amplitude ratios are independent of E and h in our approximation, the intersection of 
contours in this panel allows a direct evaluation of iV and K s . The other panels can 
then be used to determine E and ho- Therefore, observing a neutron star from the 
equator tells us more, in this context, than observing it from the pole. Note that the 
results are sensitive to the measurement of the width ratio and amplitude ratio, since 
the widths and amplitudes of the two spectra depend on the stratification and buoyancy 
in a similar manner. 



7. Comparison with Ion Collider Experiments 

The gravitational wave "experiments" proposed in §6 to measure K, N, and E dovetail 
with recent progress on measuring the compressibility and viscosity of bulk nuclear 
matter in terrestrial experiments. Indeed, it has been argued [39] that gravitational 
wave observations of pulsars will become the main experimental avenue for probing the 
nuclear equation of state in the future. 

The compressibility is usually parameterised in terms of the compression modulus 
k, which has the dimensions of energy; in our notation, we have K = Am p gL/ k (where 
A is the mean atomic number, and m p is the proton mass). It can be measured in one of 
two ways: (I) heavy-ion collision experiments, and (II) nuclear resonance experiments. 
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Figure 7. Spectral line width and amplitude contours for the plus and cross 
polarisations seen by an equatorial observer, showing the variation with N, K s and E. 
Top-left: N versus K s for E = 1CT 7 , M [h {E )+] T [fy E )+] /ho = 0.025,0.05,0.075,0.1 
(dashed contours), M [h {E)x ]T [h {E)x ] /h = 0.025,0.075,0.125,0.175,0.2 (dash- 
triple-dotted contours), T [/i(E)+] (solid contours)= T [/i(E)x] (dash-dotted contours)= 

0. 0001.0.0002.0.0005.0.001.0.002.0.005. Top-right: E versus K s for N = 

1, M [/i(E)+] r [/i(E)+] /^o (dashed contours)=M [/i(E)x] T [/i(E)x] /^o (dash-triplc- 
dotted contours)= 0.0025,0.01,0.025,0.05,0.075,0.1, T [h {E)+ ] (solid contours) = 
r[/i( E )x] (dash-dotted contours) = 0.0001,0.0003,0.001,0.003,0.01,0.05. Bottom- 
left: _E versus N for i\T s = 1, contours as for top-right. Bottom- 
right; M [ft( E )+] r [ft( E )+] /M [/1(e) x] r [ft(E)x] = 0.5,1,1.5,2,5 (dashed contours), 
T [/i(E)+] /r [ft(E) x] = !. 1-1, i- 2 , 1-3, i- 4 , 1-5 (solid contours). 

In (I), the ratio of K + meson multiplicities in Au + Au and C + C collisions (interpreted 
with the aid of nonequilibrium quantum molecular dynamics simulations) points to a soft 
equation of state, with k « 200 MeV. [66, 67, 42, 44, 43] In (II), k is deduced from the 
centroid energies of the isoscalar giant-monopole and isovector giant-dipole resonances 
in heavy nuclei. For Zr and Pb, the data yield 243 Mev < k < 268 MeV, interpreted in 
the context of a relativistic mean-field theory of the nucleus. [36, 38, 42] Again, a soft 
equation of state is preferred. In our notation, (I) and (II) imply 0.73 < K < 0.97. 

The viscosity of bulk nuclear matter is usually parametrised in terms of the shear 
viscosity to entropy density ratio, rj/s (SI units: K s). In our notation, we have 
E = (A'k B /m p L 2 il)(r]/s) for the Ekman number, where 1 < A' < 2 is the entropy per 
nucleon in units of Boltzmann's constant, fce- The viscosity to entropy ratio i]/s can be 
measured in relativistic heavy-ion collision experiments, by looking at the percentage 
azimuthal anisotropy of the hadron flux released during the collision. [68, 44] During the 
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first ~ 10~ 23 s after the collision, the nucleons flow hydrodynamically (and elliptically), 
and rj/s is proportional to the resulting diffusion coefficient, which can be related 
to the azimuthal anisotropy. Data from the latest PHENIX experiments (Au + Au 
collisions) imply rj/s ~ h/4irkB for energies just above the quark deconfinement 
phase transition. The latter value is the conjectured quantum lower bound for an 
ideal fluid. It has been derived by exploiting the duality between a large class 
of gauge theories similar (but not identical) to quantum chromodynamics and the 
thermodynamics of black hole horizons. [69, 45] In our notation, the measurements 
imply E = 5 x \Q- l7 A\L/W m)- 2 (fi/rads" 1 )- 1 . Interestingly, this is ~ 10 orders of 
magnitude smaller than the Newtonian viscosity for neutron-neutron scattering in a 
superfluid calculated by Cutler and Lindblom. [70, 56] 

8. Conclusions 

In this paper, the gravitational wave signal emitted by the nonaxisymmetric spin-up 
flow excited by a pulsar glitch is calculated analytically in the context of a cylindrical, 
Newtonian-fluid toy model, which neglects the physics of pinning and turbulence in 
the boundary layer between the superfluid and the crust. The model clarifies the 
physical roles played by stratification and compressibility in nonaxisymmetric spin up. 
It is found that stratification and compressibility compete to increase and reduce the 
buoyancy respectively, and consequently, reduce and increase the spin-up time. Higher 
order mode numbers enhance the effects of buoyancy. Equations (68), (71) and (72) 
give approximate formulas for the gravitational wave spectrum in the plus and cross 
polarizations. Measurements of the amplitudes and widths of the lines in these spectra 
can be used to constrain K S) N, and E, by reading off the intersection of the contours 
plotted in Figures 6 and 7. 

Current estimates of the internal properties of a typical neutron star imply that 
the characteristic wave strain (64) produced by a large (SQ/Q > 10~ 4 ) glitch evaluates 
roughly to h > 5 x 10~ 26 for L = 10 4 m, mass M* = M , g « GM^/L 2 , d = 1 kpc, and 
/ = Q/2tc > 0.2 kHz. While this estimate is crude, in view of the over-simplified model 
assumptions discussed in §1 and §2, it suggests that a glitch with these parameters 
arguably approaches the threshold of detectability by a coherent, matched filter search 
(1% false alarm, 10% false dismissal) with Advanced LIGO for an observation time of 
T obs = 14 days. Note that T obs is limited by two factors: (1) the computational expense 
of a fully coherent search with matched filters, which is prohibitive for T Q b s > 14 d at 
present [71]; and (2) the duration of the exponential relaxation phase following a glitch 
(i.e., the Ekman time-scale), which is observed to last from hours to a few weeks in radio 
pulsar timing experiments [13, 46, 47]. Restriction (2) is lifted if the nonaxisymmetry 
persists even after the relaxation phase, e.g. due to inhomogeneous vortex pinning 
[32, 64]. Of course, compressibility and stratification act to reduce h+ and h x below the 
fiducial estimate of h given above, sometimes substantially; one has 1 < N < 10 and 
0.1 < K s < 1 for a typical neutron star [17, 35, 70, 56]. The above estimates apply to 
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the largest glitches observed to date, with 5Q/Q > 10~ 4 [11]. A recent analysis of the 
observed glitch size and waiting-time distributions in individual pulsars suggests that 
future monitoring will reveal even larger glitches [12]. Once detected, the glitch signal 
will allow us to compare astronomical measurements of K and E with data from recent 
heavy-ion collider experiments. The collider data prefer a soft nuclear equation of state 
and a viscosity close to the quantum lower bound for an ideal fluid, contradicting some 
neutron star models. 

Although our toy analytical model involves many unavoidable approximations (e.g. 
cylindrical geometry, uniform gravity, Newtonian fluid, idealized superfluid boundary 
layer), it nevertheless demonstrates in principle for how the gravitational wave signal 
from a pulsar glitch can be inverted to extract important information about neutron 
stars and bulk nuclear matter, independently from terrestrial experiments with particle 
accelerators. Our preliminary detectability estimates offer grounds for hope that such 
gravitational wave observations are realistic in the medium term. Of course, to apply 
these ideas in practice, the idealised model in this paper is inadequate and must be 
supplanted by large-scale numerical simulations, e.g. Refs [18, 23]. 
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